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ABSTRACT 



Two stationary first-order autoregressive processes with Beta marginal 
distributions are presented. They are both linear, additive processes but the 
coefficients are Beta random variables. Their autocorrelation functions are 
investigated: one is positive and the other alternates in sign. The useful- 

ness of the models in simulation is discussed. The Bivariate Beta distribu- 
tions of two consecutive observations are considered in some detail. Several 
examples are given, including a Bivariate Uniform process which is also 
examined in detail. The relationship of these Bivariate Beta distributions 
to the Dirichelet distribution is discussed. 
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1. INTRODUCTION 

The Beta distribution is the most versatile and useful distribution 
available on a bounded interval. Despite this there are very few practical 
models for describing correlation between pairs of Beta random variables or 
serial correlation in sequences of them. This is unfortunate since there is a 
natural interest in modelling sequences of dependent Beta variates. These 
arise in a variety of ways but are often associated with the stochastic 
behaviour of a proportion or probability over time. One example is the study 
of market share, e.g. the proportion of the market held by a particular 
product. (See, for example, Wichern and Jones (1977)). 

Some efforts have been made to model such behaviour. Azzalini (1982) 
developed a simple Markov model for use in a quality control context where the 
proportions defective in adjacent batches sampled may be expected to be depen- 
dent. He used the Product Autoregressive process discussed by McKenzie (1983). 
This model has limited usefulness for the Beta distribution, however, since it 
requires that of the two parameters of that distribution one must be integral 
and the other rational. Another approach has been suggested by Souza and 
Harrison (1981). They attribute to the proportion a Beta distribution which is 
revised with each observation by a procedure based on Bayesian and Information 
Theoretic concepts. This approach appears to offer some promise for forecast- 
ing but is fairly complex in general. A more traditional time-series approach 
to such a data set would be to treat the proportions not as Beta's but as 

Logistic-normal variates. Thus, we would transform to log-odds, i.e. 

Y = Jln[X /(1-X )] , and attempt to model {Y } as an autoregressive moving- 

L L L U 

average process. This approach is suggested by the work of Aitchison and Shen 
(1980) and Aitchison (1982), but does not appear to have been investigated as 
a time-series procedure yet. It discards the Beta marginal distribution which 
is our main interest here. Further, it appears to be extremely difficult to 
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reverse the procedure and generate a sequence of proportions with a specific 
correlation structure. 

The purpose of this work is to present a discrete time Markov process with 
a Beta marginal distribution. It is constructed in the spirit of the recent 
work on modelling of non-6aussian time series illustrated, for example, by 
Lawrance and Lewis (1980, 1982) and Jacobs and Lewis (1983). Since simulation 
is also a major motivation for such work, we seek models which are simple and 
flexible and whose parameters are few and physically meaningful. The aspect of 
simulation is important here for, as Schmeiser and Lai (1980) noted in a recent 
survey (1980), there are few practical ways of generating dependent Beta random 
variables. This is because the usual multivariate Beta distributions are 
closely related to the Dirichelet distribution and so constrain their vector 
variates in a way which is undesirable in general. 

We present here two simple discrete-time stochastic processes whose mar- 
ginal distributions are Beta random variables. The processes are linear, 
additive autoregressions with random coefficients. The coefficients themselves 
are also Beta variates. There is a single free parameter which corresponds in 
a simple way to the correlation in each model and a wide range of correlation 
is possible. Because of their simplicity, the models provide a powerful way of 
generating sequences of dependent Beta variates using only independent Betas. 

As noted above, there is a scarcity of practical multivariate Beta distribu- 
tions. Thus, the bivariate distributions associated with the processes are 
discussed in some detail. They exhibit a number of interesting features. We 
also examine in detail the particular case of the bivariate Uniform distribu- 
tion. It plays an important role in the simulation of dependent pairs of 
random variables. 
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Before describing the models we may not that there is a well-known 
continuous-time Markov process with Beta marginals. It arises in genetics and 
is one of the forms of the Wright-Fi sher gene frequency models described in 
detail in the books by Karlin and Taylor (1975, 1981). It is a diffusion 
process and, in different forms, has found applications in sociology, psychol- 
ogy and marketing. It is also derived by Massey et al. (1970) as a stochastic 
response model. They obtain it as the limiting form of the "contagious 
binomial" distribution developed by Coleman (1964) to model voting behaviour. 
Of course, such processes are in continuous time and it is by no means clear 
how they can be restructured in discrete time. Nor is it clear whether the 
processes presented here represent some discrete time formulations of the 
diffusion processes. 

2. THE MODELS 

2.1. A random variable (r.v.) X is said to have a Beta distribution with 
parameters (a,B) if it has probability density function (p.d.f.) 

x*^~^( 1-x)^”^ 

" B(a,e) , 0 < X < 1 ; a,B > 0 . 

For convenience in what follows, we shall write such a random variable as 
Be(a,6) . We note for later use that for such X , 

E(X) = a/(a+g), Var(X) = ag/C(ot+g) (a+g+1)] and the third moment about the 

3 

mean is given by m^ = 2ag( 6-a)/[(a+g) ( a+g+1) ( a+g+2) ] . 

The models presented here use the following results; 

1 - Be(a,g) = Be(g,a) (1) 



Be(a,g) • Be(a+B,y) = Be(a,6+T) 



( 2 ) 
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The first of these two results is well known and easily demonstrated. The 
second result states that the product of two independent Beta r.v.s. with 
parameters as specified is itself a Beta r.v. The result may be verified by 
considering the Mellin Transform (Widder, 1946), i.e. E(X^) , of the product 
on the left-hand side of (2). It is 

p/wSv _ B(g+S ,B) ^ B( g+g+S ,y) _ T( crj~s) T( g+g+y) _ B( g+S , g+y) 

' B(g,g) B(g+g,y) T( g) T( g+g+y+S ) B(g,g+y) 

which is the transform of the right-hand side of (2). We shall refer to the 
application of (2) to change one Beta r.v. into another as the Beta-Beta 
transformation. 

Two distinct models are presented here. One is for positively correlated 
pairs of Beta r.v.s. and is denoted by PBAR and the other is for negatively 
correlated pair's of Beta's and is denoted by NBAR. Both models are linear and 
additive and have random coefficients. 

2.2. The PBAR model is given by 

Xt = 1 - U^d - Vt-1> (3) 

where {U^} and {W^} are independent sequences of independent identically 
distributed (i.i.d.) r.v.s., independent of previous X’s , and is 

Be(g,g-p) and is Be(p,g-p), (0 < p < g) . Now if X^_^ is Be(g,g) then 

\^t-l Be(p,g+g-p) from (2) and 1 - W^X^_^ is Be(g+g-p,p) from (1). 

Further use of (2) shows that U^(1-W^X^ is Be(g,g) and so X^ given by (3) 
is Be(g,g) by (1). Thus, equation (3) defines a stationary process {X^} with 
a Beta(g,g) marginal distribution. 

Further, there is a single free parameter in this scheme, viz. p . As we 
shall see, the value of p determines the correlation structure of the 
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process . From the structure of (3), the process is a first-order linear 

autoregression with random coefficients. Direct calculation yields the auto- 
correlation function of the process as Pj^(k) = , k = 0 , 1 ,..., where 

p = E(U)E(W) = pe/a(a+6-p) . (4) 

Now p as defined by (4) is a monotonic increasing function of p for fixed 
(a,e) . Further, since e>0 ,0<p<a,we may deduce that 0 < p < 1 . 

Thus, the entire range of positive correlation is possible for any values of a 
and e , i.e. any Beta marginal distribution. 

There are two limiting cases for the parameter p we may consider. When 
p is zero = 0 with probability one and is Be(6,a) . Thus, is 

independent of ^ and p = 0 . When p = a both and W^ are unity 

with probability one. Thus, 1 ^nd p = 1 . The process is not ergodic 

in this case. 

It is important to notice that with this model (3) p = 0 implies 
i ndependence . 



2.3. The model for the NBAR process is given by 



X. = 



= v^d - 



Vt-1> 



(5) 



As before, {V^} and {W^} are independent sequences of i.i.d., r.v.'s 
independent of X^ ^ and is Be(a,6-p) and W^ is Be(p,a-p) . 

Again, it is easily verified using (1) and (2) that equation (5) will generate 
a stationary process whose marginal distribution is Beta(a,g) . The NBAR 
process is also a linear autoregression with random coefficients. It too has 

k 

autocorrelation function of the form p , k = 0,1,..., but now 



p = - p/ ( a+B-p) • 



( 6 ) 
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Notice that the specification of the distributions of and requires 

that 0 < p < a and 0 < p < e . From (6), p is a monotonic decreasing 
function of p for fixed (a, 3), and so we may deduce that for the NBAR 
process{5) -max(a/3,6/a) < p < 0 . 

The upper extreme of zero is again attainable when p = 0 and > 

which is independent of ^ . As with the PBAR process, it is important to 
note that p = 0 implies idependence. The lower limit -max{a/S,3/a) is also 
attainable. If 3 < a it is attained when p = 3 and so = 1 . If 

3 > a it occurs where p = a which corresponds to = 1 . When a = 3 

the lower limit is -1 which corresponds to the usual antithetic relationship, 

= 1 - , given by = 1 . 

2.4. The models PBAR and NBAR given by (3) and (5) yield random coefficient 
autoregressions of order one. Further, the first-order correlation p satisfies 

- max(-| , -^) < p < 1 . (7) 

This is not the greatest possible range for general (a, 3) . For example, if 
a = 2, 3 = 1 we find that -0.5 < p < 1 for these models. On the other hand, 
it may be deduced from Moran (1967) that the correlation between two Be(2,l) 
r.v.s. is bounded below by (9 t:- 32)/4, i .e. approximately -0.9314. How far the 
lower bound given by (7) is from the minimum correlation possible is not known for 
general (a, 3) . 

We may note, however, that in the symmetric case, i.e. a = 3 , the range 
is (-1,1) as would be hoped. As before, the two extremes may be attained: 
the upper limit of 1 from ~ ^t 1 lower limit of -1 from the 

usual antithetic relationship ^ 1 * latter is obtained from 

NBAR with p = a = 3 SO that = 1 . In particular, it is now 
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possible to generate r.v.s. uniform on (0,1) with any first-order autoregressive 
correlation, i.e. any pe[-l,l] . We return to this point later. 



3. BIVARIATE DISTRIBUTIONS OF (X^, X^_^) 

3.1. For convenience, we rewrite the PBAR model (3) as Y = 1 - U(l-WX) . 
Initially, consider the joint p.d.f. of (Y,X) conditional on U . It can be 
written in the form 






t^(x) 



f 

W^ XU 



1 - y < u < 



ilZ 

1-x 



The Joint p.d.f. of Y, X and U may now be derived and from it the p.d.f. 
of (Y,X) is obtained in the form 



T ,A Q 



( 8 ) 



where m(y,x) = min(y ,x(l-y)/(l-x) ) and = B(a,B)B(p,a-p)B(B,a-p) . 

A change of variable t = (l-x)s/(l-y) in the integral in (8) yields the 
result fy j^(y,x) = fy j^(x,y) . Thus, the bivariate p.d.f., which is defined 
on the unit square, is symmetric about y = x . Note also that m(y,x) = y 
if y < X . Thus, the two forms of the integral exist on either side of 
y = X , and we can define 

g(y,x) y < X 

g(x,y) y > X 

where g(y,x) is given by the right hand side of (8) with m(y,x) replaced 

by y . 
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This symmetry is also important from the viewpoint of modelling or iden- 
tification of the PBAR process. The importance arises from the idea of time- 
reversibility of a stationary process. The concept is discussed in detail by 
Weiss (1975). A discrete-time stationary process {X^> is time-reversible if 
the joint distributions of {X^ . ,X^} and {X^,X^_^,...,X^} are identical 

for every t . In the case of a first-order autoregression , as here, the 
joint p.d.f. can be expressed as 

t-1 t-1 

f(x, ,x„, ... ,x ) = n f(x.,x. ,)/ n f(x.) 

1 ^ i = l 1 i = i 1 



Such a process is time reversible then whenever the bivariate distribution is 
symmetric. Thus, the PBAR process is time-reversible. 

Indeed, the time reversed process is given by 



X 



t 






) 



where (U'^} and {Wp have exactly the same properties as {U^} and {W^} 
respectively defined by (3). 



3.2. For the NBAR process, writing Y = V(l-WX) and proceeding as for the PBAR 
process yields 






/, \B-1 m(l-y,x) 1 1 „ 1 

(1-x) f ' p-l|- ,, > no-p-l/i \S-p-lj 

= J s*^ [xy - (l-x)s] ^ (1-y-s) ^ ds 



( 10 ) 



where C_ = B(a,B)B(p,a-p)B(a,S~p) . 

Again, the change of variable t = (l-x)s/(l-y) shows that the bivariate 
p.d.f. is symmetric about y = x and the NBAR process is time-reversible. 

The structure of the density, however, is a little more complex than for the 
PBAR process. The upper limit of the integral in (10) is m(l-y, x) = 1 - y 



9 

if X + y > 1 . Thus, the form of f depends upon which side of the line 
X + y = 1 we are on. 

Define 

9i(y»x) X + y < 1 , 

’ 92(y»x) X + y > 1 . 



From symmetry about y = x , we know that g^(y,x) = g^(x,y) and 
92 (y»x) = 92 (^»y) • 0^ immediate interest is the relationship between 

g^ and g^ . A further change of variable in the integral in (10) yields the 
following result. Using an obvious notation to denote dependence upon the 
parameters of the marginal distribution: 



g^(l-y, 1-x; a, 6) = g 2 (y, x; g, a) (12) 

Using (12) to evaluate f~ in both the triangular regions induced by the line 
X + y = 1 yields 

*x,0) *’Ty^(y,x, g,ci) (13) 



This result (13) is an obvious two-dimensional analogue of the relationship 
specified by equation (1), viz. f^(l-x; a,g) = * 

Further, by the symmetry about y = x , equation (12) yields 

g^(l-x, 1-y; a,g) = 92 (y.x; g,a) (14) 

which specifies the nature of the relationship between g^ and g^ across 
the line x + y = 1 . In particular, note that when a = g the bivariate 
p.d.f. is symmetric about both y = x and x + y = 1 . 
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4. MOMENTS 

Although it is difficult to obtain the p.d.f.s in an explicit form, much 
information can be derived directly from the structural relationships, (3) and 
(5). This is also true of the conditional distributions. In both cases, the 
conditional expectation is linear. For the PBAR process, using (3) yields 

E(Xt^l|Xt = x) = [(a+e-p-ag) + pgx]/a( o+0-p) , (15) 

and for the NBAR process, (5) yields 

E(X^.^l|Xt = x) = (a-px)/(a+g-p) . (15a) 

Further, the time-reversibi 1 ity of the processes ensures that the inverse 

regressions are identical to (15) and (15a), i.e. E(X^lX^^^ = x) = E(X^^^|X^ = x) 

In both cases the conditional variances are quadratic and of the form 

2 2 2 2 2 2 2 

var(X^^^|X^ = x) = a (1-njx) + (cr +m )a^jX where = var(W) and m and 

2 

a are the mean and variance of U for PBAR and of V for NBAR. 

Higher order moments are important in the identification of non-standard 

time-series models and we note here two of particular interest. The first is 
'><2 ^ 

C 2 j(k) = Cov(X^, For time reversible processes C 2 ^(k) = . 

For both the PBAR and NBAR processes, C 2 ^(k) = ^ ~ 0,1,2,..., where 

m^ is the third moment of X^ about its mean, y say, and these are given 
in Section 2.1. 

Another moment of particular value in residual analysis is based on the 

residuals from minimum mean square error prediction, i.e. 

R^ = (X^ - u) - - u) . Such residuals are uncorrelated for PBAR and 

2 

NBAR processes. Of more interest is the behaviour of Cov(R 2 , R^ |^) which 
is useful in distinguishing between constant and random coefficient models. 

For both the PBAR and NBAR processes 
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m3(l-p)(l-p^)p^ k = 1,2,... 

Cov(R^, R^_^) =0 k = -1,-2,... 

m3(l+2p)(l-p)^ k = 0 . 

Further details and discussion of the usefulness of these moments may be found 
in Lewis and Lawrance (1983b). 

5. EXAMPLES 

To illustrate the nature of the bivariate p.d.f. the functions g from 
(9) and and g^ from (11) are evaluated explicitly for a few specific 

cases below. 

PBAR 

PI: a = B = 2; p = 1 i .e. p = 1/3, 

g(y,x) = 12y(l-x) 

P2: a = 3, B = 2; p = 1 i .e. p = 1/6, 

g(y,x) = 72(l-x)C(l-x)y(2-y)-2(l-y)tn(l-y)-2y(2-x-y)] 

P3: a=4, b=2;p= 3 i .e. p = 1/2, 

g(y,x) = 120(l-x)[l+2(l-y)in(l-y) - (1-y)^] 

NBAR 

Nl: a = B = 2; p = 1 i.e. p = - 1/3, 

g^(y,x) = 12xy 
g^(y,x) = 12(l-x)(l-y) 

N2: a = 3 , B = 2; p = 1 i .e. p = - 1/4 

g^(y,x) = 36x^y^ 

g^(y,x) = 36(l-x)(l-y)(x+y + xy-1) 

N3: a = 5, B = 13/3; p = 4 i.e. p = - 3/4 

g 2 (y»x) = A(l-x) (1-y) 

9n(y,x) = g 2 (y,x) - B(l-x-y)^/^3(l-y)^(l-x)^ _ |( l-y)2(i.x)2(i_x.y) 

+ y(i-y)(i-x)(i-x-y)^ - 
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where A = 3532100/2187 = 1615.04, and B = 140A/283 = 798.96. 

Plots of the contours and the surface of PI are displayed in Figures 1(a) and 

(b). If they are rotated through 90° about the point (0.5, 0.5) the corre- 
sponding plots are obtained for N1 . Contour plots are displayed for P2 and N2 

in Figures 2 and 3. 

Figures 1,2,3 

The most remarkable feature of these distributions is the appearance of a 
ridge in the surface of PI (and so Nl) but not in P2 or N2. The ridge is due 
to the fact that both densities have two forms as described by equations (9) 
and (11). Thus the ridge, if it occurs, corresponds to the line y = x for 
PBAR p.d.f.s and x + y = 1 for NBAR's. From the definition, the p.d.f.s will 
be continuous but tneir derivatives need not be so. Any points of disconti- 
nuity occur on these two lines. To determine the general conditions for 
occurrence of a ridge we can examine the behaviour of the derivatives of the 
p.d.f.s on both sides of the two lines. This procedure yields the following 
results. For the PBAR density (8) there is no ridge provided p < (a-1) . 

Using (4), we may deduce that no ridge occurs provided p satisfies 

0 p < B( a-l)/a( B+1) . (16) 

For the NBAR density (10) there is no ridge provided p < (a+B-2)/2 i.e. 
provided 



0 ^ p > [2 - (a+g) ]/(2+crt-0) . (17) 

These conditions (16) and (17) are violated by P3 and N3 respectively and so 
both densities have a ridge. These are illustrated in the surface plots of 
P3 and N3 displayed in Figures 4 and 5, respectively . 



Figures 4,5 
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6. THE UNIFORM PROCESS 

A bivariate distribution of particular interest is that with Uniform 
marginals. Apart from the natural interest in modelling data from such a dis- 
tribution it has an important application in simulation. By using the inverse 
distribution function transformation of the Uniforms bivariate distributions 
with any other marginals can be obtained. This is an important approach to 
the generation of pairs of dependent random variables. 

Some recent development of Bivariate Uniform distributions appears in the 
papers by Barnett (1980) and Lewis and Lawrance (1983a). The former exhibits 
several different Bivariate Uniform distributions but they are generally 
complex, have limited correlation, and require distribution function trans- 
formations to obtain the Uniforms. The latter work gives several procedures 
for generating a pair of dependent Uniforms from a random coefficient regres- 
sion on independent Uniforms. The procedures are simple and the entire range 
of correlation can be attained. However, to achieve this breadth the correla- 
tions are usually complex functions of two parameters. 

The Bivariate Uniform with p.d.f.s given by (8) and (10) is particularly 
useful in this kind of application. The entire range of correlation [-1, 1] is 
available to the Uniforms and so the entire possible range will be available 
to the transformed variates. Further, the correlation is a simple function of 
a single parameter, i.e. p = ± p/(2-p) , (0 _< p <_ 1) and so any desired cor- 
relation is easily achieved. Finally, the generation of the desired pair is 
straightforward involving only the additional generation of two independent 
Beta r.v.s. 

Since the Bivariate Uniform corresponds to (8) and (10) with a = B = 1 
considerable simplification is possible. In this case, 

Ty ~ x(l-y»x) 3nd so only f^ is considered here. Making a suit- 

able change of variable in the integral yields a more useful expression for 
the density given by (9), viz. 
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g(y,x) = 



2(1-P) 



B(p,l-p)xP(l-y)P 0 



(tane) 



2p-l 



[1 



y(i-x) 
■ x(l-y) 



2 

sin e] 



-P 



de, (y < x) 



The density has a line singularity on y = x as can be seen from the contour 
and surface plots shown in Figure 6 for the case p = 1/2 . 



Figure 6 



7. RELATIONSHIP WITH THE DIRICHELET DISTRIBUTION 

The best-known "Bivariate Beta" distribution is the Dirichelet distrib- 
ution. It is described in detail by Johnson and Kotz (1972). It is defined 
not on the Unit square as are the densities given by (8) and (10) but on the 
Unit simplex, i.e. {(x,y): 0<x, y<l,x+y<l>. As such it plays a 

natural role as the joint distribution of two proportions from a single popu- 
lation. It seems to be generally regarded as a Bivariate Beta because both 
the marginal and conditional distributions are Beta. The joint p.d.f. in the 
case of identical Be(a,B) marginals is 



f(y»x) = 



r(g+B) 

r(a)r(a)r(B-a) 



a- 1 a- 
X y 



^(i-y-x)^' 



x-l 



and the correlation between X and Y is p = - a/6 (a < b) . Note that p 
here depends explicitly upon the parameters of marginal distribution and so is 
fixed for any particular marginal distribution. Further, -a/6 is the minimum 
correlation attainable from the NBAR process. As noted, it is attained when 
p = a and Y = V(l-X) where V is Be(a,6-a) and independent of X . 

Thus, Y/(l-X) is independent of X , and symmetry ensures that X/(l-Y) is 
independent of Y . This characterizes (X,Y) as having a Dirichelet distrib- 
ution by a result of Darroch and Ratcliff (1971). Thus, we may view the 
Dirichelet distribution as a limiting form of the NBAR process distribution, 
as p ^ min (a, 6) . 
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8, EXTENSIONS 

It is possible to extend the time-series models (3) and (5) to higher 
orders of dependence. However, two simple and more immediate extensions lie 
closer to the area of simulation and are noted here. The fact that the PBAR 
and NBAR models yield simple but powerful methods of generating sequences of 
correlated Beta r.v.s. has been emphasized. Such sequences are stationary so 
that each Beta r.v. has the same distribution. An obvious extension is to the 
generation of pairs of dependent Beta variates with different distribution. 
This may be achieved in a variety of ways using (1) and (2) and (3) or (5) if 
we wish. 

A second simple generalization is to bounded intervals other than (U,l). 
Since an alternative sample space is achieved by a linear transformation the 
procedure is straightforward and all correlations are unaffected. Thus, sup- 
pose we wish to develop a PBAR process for Be(a,B) r.v.s. defined on (a,b) 
rather than (0,1). By considering the X^'s transformed to Y^'s on (0,1) 
and {Y } satisfying (3) we find that the PBAR for {X } is given by 

X^ = b - U^[b-a-w^(x^-a)3 . 
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Figure 1(a). Bivariate Beta PI, p=0.33, contours ,5(.5)2.5 



17 




Figure 1(b). Bivoriate Beta PI, p=0.33, Surface Plot 



Figure 2. Bivariate Beta P2, p=0.1667, contours .5(.5)3 
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Figure 3. Bivariate Beta N2, p=-0.25, contours .5(.5)3 
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Figure 4. Bivariate Beta P3, p=0.5 
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Figure 5. Bivariate Beta N3, p = -0.75 



Figure 6. Bivariate Uniform, p=0.333, contours .5(.5)3 
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Figure 6. Bivariate Uniform. p=0.333, contours .5(.5)3 
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